 function [logPDFNoConstant,logPDF]=MNIWPDF(X,W,muMat,Pprec,S,v)
% 
%  [logPDFNoConstant,logPDF]=MNIWPDF(x,W,muMat,Pprec,v)
% 
% Computes the Log pdf of (X,W)~MNIW(muMat,P,S,v)
% such that 
% X|W ~ MN( muMat, W   kron P ) 
%   W ~ IW( v    , S          ) 
% 
% See matricvariateNPDF.m and invWishartPDF.m 
%     for the computation of the invidual densities 
% 
% *Note* uses inv(P) instead of P 
% 
% Notation follows from Lubrano et. al page 308.
%
% The unormalized pdf will not be computed --hence no need to account for constant--
% if only one output argument is requested. 
%
% Code uses a SVD to obtain the inverse and determinant of W 
% 
%% Input 
% X:     [p,q] matrix
% muMat: [p,q] matrix with Mean 
% Pprec: [p,p] matrix with *Prec=INV(P)*
% S:     [q,q] matrix for IW 
% v:     (scalar) degrees of Freedom for IW 
%
%% Output 
% logPDFNoConstant: log PDF without the normalizing constant 
% logPDF          : log PDF with normalizing constant 
% 
% Alejandro Justiniano  Aug 30 2013 (C) 
% =========================================================================
[Nr,Nc]=size(X); 
% Nr=k 
% Nc=n 
%% 1. Obtain inv(W) using the SVD and use this to compute the log determinant 
W=0.5*(W+W'); 

% This is more robust but probably slower than inv(W)
% for n small 

%% 1.a SVD 
% PP*DD*PPinv'=W  Notice the transpose 
% PPinv=inv(PP)' 
% PPinv'=inv(PP); 
[PP,DD,PPinv]=svd(W); 

%% 1.b Truncate small singular values 
tolZero=1e-12; 
firstZero=min(find(diag(DD) < tolZero));
if isempty(firstZero)==false 
    PP=PP(:,1:firstZero-1); 
    PPinv=PPinv(:,1:firstZero-1); 
    DD=DD(1:firstZero-1,1:firstZero-1); 
end 

%% 1.c Inverse and determinant 
Winv=PP*(DD\eye(Nc))*PPinv'; 
logDetW=sum(log(diag(DD))); 

%% Check inverses and determinants are correct Aug 31 2013 
% fprintf('Max diff in Winv %1.12f \n', ...
%     max(max(abs( Winv -inv(W) ) )) );
% fprintf('Max diff in log(det(W))  %1.12f \n', ...
%     abs(logDetW-log(det(W))) );

%% 2. Log PDF w/o Normalizing constant
% Note the change from + to - in the second line 
% since using the log(det(W)) not log(det(W^-1)) 
% see matricVariateNPDF.m
logPDFNoConstant=-0.5*trace( Winv*(X-muMat)'*Pprec*(X-muMat) )-...
    0.5*Nr*logDetW-0.5*trace(Winv*S)-0.5*(Nc+v+1)*logDetW; 

if nargout < 2 
    return 
end 
%% 3. Log Constant 
logConstantIW=zeros(4,1); 
logConstantIW(1)=v*Nc*log(2)/2; 
logConstantIW(2)=Nc*(Nc-1)*log(pi)/4; 
logConstantIW(3)=sum(gammaln(((v+1 -(1:Nc)))/2)); 
logConstantIW(4)=(-v/2)*log(det(S));
logConstant=-sum(logConstantIW)-0.5*Nr*Nc*log(2*pi)+...
    0.5*Nc*log(det(Pprec)); 

%% 4. Log PDF with constant 
logPDF=logPDFNoConstant+logConstant;

%% Check Densities without and with constants match 
%  Aug 31 2013 
% checkNC=zeros(2,1); 
% checkCONS=zeros(2,1); 
% [checkNC(1),checkCONS(1)]=matricvariateNPDF(X,muMat,inv(W),Pprec); 
% [checkNC(2),checkCONS(2)]=invWishartPDF(W,v,S); 
% if nargout < 2 
%     return 
% end 
% fprintf('Max diff in logPDFMNIWX  w/o constant %1.12f \n', ...
%     abs(logPDFNoConstant-sum(checkNC)) );
% fprintf('Max diff in logPDFMNIWX  w  constant %1.12f \n', ...
%     abs(logPDF-sum(checkCONS
